Importação da base de dados

dataset_diario é a base com valores de fluxos diários de 2002 a 2009.
dataset_mensal é a base com os valores de fluxo mensais de 1992 a 2009.

Nota: a informação sobre os dias da semana é redundante, pois pode ser obtida através de:

dia <- weekdays(index(dataset_diario))

Por isso, na leitura da base em dataset_diario, esta coluna foi ignorada.

Séries Temporais

Conversão dos dados para o formato de série temporal no R.

Série Mensal:

# Frequency --> número de observações por unidade de tempo 
# define a unidade de tempo (e.g. 12: unidade de tempo = ano)
tsMensal <- ts(dataset_mensal, frequency=12, start=1992)

plot of chunk unnamed-chunk-6

Série Diária:

Série temporal diária

tsDiario <- msts(dataset_diario, start=c(2002,1,1), seasonal.periods=c(7, 365.25))

plot of chunk unnamed-chunk-8

Nota: é importante observar os anos bissextos, 2004 e 2008.

Inspeção dos Dados

Série para um dia específico da semana

$plot of chunk unnamed-chunk-9$

Série para um dia e mês específico da semana

plot of chunk unnamed-chunk-10

Observar a influência dos feriados. Exemplo, 12 de outubro foi uma quinta-feira em 2006.

Inspeção dos Dados - 2

A série temporal pode ser decomposta em:

Trend: padrão de crescimento ao longo do tempo.
Seasonal: padrão recorrente em um determinado período de tempo.
Irregular (random): parte irregular (ruidosa), sem padrão sazonal ou tendência.

Fluxo Mensal

plot of chunk unnamed-chunk-11

Fluxo Diário

plot of chunk unnamed-chunk-12

Observações

Clara tendência (trend) de subida.
Sazonalidade de 1 ano.
No fluxo diário, há sazonalidade semanal.

Inspeção dos Dados - 3

Fluxo Diário para um intervalo arbitrário de 8 semanas

$plot of chunk unnamed-chunk-13$

Inspeção dos Dados - 4

Auto-correlação

acf(dataset_diario, col="red")

plot of chunk unnamed-chunk-14

Parcial:

pacf(dataset_diario, col="red")

plot of chunk unnamed-chunk-15

Objetivos

Realizar prospecções de curto, médio e longo prazo.

Curto Prazo

30 Dias
45 Dias

Para curto prazo, será utilizada a série temporal com dados diários (tsDiario).

Médio Prazo

4 Meses
6 Meses

Longo Prazo

1 Ano
2 Anos

Para médio e longo prazo, será utilizada a série temporal com dados mensais (tsMensal).

Separação dos dados

Separação das séries temporais em conjuntos de treinamento e de testes.

Curto Prazo (`tsDiario`)

Conjunto de Treinamento

tsDiarioTrain <- window(tsDiario, end=c(2008,183)) # até 30/06/2008

plot of chunk unnamed-chunk-17

Conjuntos de Teste

30 dias:

# De 01/07/2008 a 30/07/2008
tsDiarioTest30Days <- window(tsDiario, start=c(2008,184), end=c(2008,213))

45 dias:

# 45 dias a partir De 1/1/2008
tsDiarioTest45Days <- window(tsDiario, start=c(2008,184), end=c(2008,228))

Médio/Longo Prazo (`tsMensal`)

Conjunto de Treinamento

tsMensalTrain <- window(tsMensal, end=c(2007, 12)) # De Jan 1992 a Dez 2007

plot of chunk unnamed-chunk-24

Conjuntos de Teste

4 meses:

tsMensalTest4mth <- window(tsMensal, start=2008, end=c(2008, 4)) # De Jan 2006 a Mar 2006

6 meses:

tsMensalTest6mth <- window(tsMensal, start=2008, end=c(2008, 6)) # De Jan 2006 a Jun 2006

1 ano:

tsMensalTest1yr <- window(tsMensal, start=2008, end=c(2008, 12)) # De Jan 2006 a Jan 2007

2 anos:

tsMensalTest2yr <- window(tsMensal, start=2008, end=c(2009, 12)) # De Jan 2006 a Jan 2008

Bibliotecas Utilizada

Prospecções Estatísticas

ver (R. J. Hyndman 2014) e (Leek 2014):

library(forecast)

Redes Neurais:

library(nnet)

Estratégia de Avaliação dos Métodos

Apresenta-se a seguir:

a predição proporcionada por diferentes métodos (curva vermelha).
intervalos de predição correspondentes aos níveis de confiança 80% (em azul) e 95% (em cinza).
acurácias in-sample e out-of-sample.

Nota:

in-sample: medida no conjunto de treinamento.
out-of-sample: medida no conjunto de testes.

Acurácia

MAPE:

Diferentemente do ME, MAE e RMSE, o medidor MAPE é independente da escala.

\[M = \frac{1}{n} \sum \limits_{t=1}^{n} \left| \frac{A_t - F_t}{A_t} \right|,\]

onde $A_t$ é o valor verdadeiro, $F_t$ é o valor predito e $n$ é o número de amostras temporais preditas.

Como desvantagem, o critério MAPE tem comportamento indefinido quando o valor verdadeiro $A_t=0$ (divisão por zero).

Benchmark Métodos Estatísticos

Curto prazo

30 dias

##           Training set Test set
## Mean            14.656   22.506
## Näive            6.214    4.674
## N-Sazonal        8.720    9.066
## Drift            6.216    4.788
## Regressao        5.868    7.383

Médio prazo (4 meses)

##           Training set Test set
## Mean             25.58   35.295
## Näive             4.44    5.639
## N-Sazonal         7.68    6.051
## Drift             4.47    6.542

Longo prazo

1 ano

##           Training set Test set
## Mean            25.582   39.867
## Näive            4.440    5.708
## N-Sazonal        7.680    6.898
## Drift            4.470    4.722
## Regressao        4.741   10.023
## ETS              1.442    2.603

2 anos

##           Training set Test set
## Mean            25.582   41.232
## Näive            4.440    6.801
## N-Sazonal        7.680    9.046
## Drift            4.470    5.210
## Regressao        4.741   10.398
## ETS              1.442    2.964

Modelo auto-regressivo (AR)

Em um modelo de auto-regressão, prevemos a variável de interesse usando uma combinação linear dos valores passados da variável.

Assim, um modelo auto-regressivo de ordem $p$, denotado por $AR(p)$, pode ser escrito como:

$y_t = c + \varphi_1 y_{t-1} + \varphi_2 y_{t-2} + \dots + \varphi_p y_{t-p} + e_t,$

onde $c$ é uma constante, $e_t$ o ruído branco, e, os valores defasados de $y_t$ são os preditores.

Para um modelo $AR(1)$, tem-se $\varphi_1\in(-1,1)$

Para um modelo $AR(2)$, tem-se $\varphi_1,\varphi_2\in(-1,1)$ com $\varphi_1 + \varphi_2 < 1$ e $\varphi_2 - \varphi_1 < 1$

Modelos médias-móveis (MA)

Ao invés de usar valores passados da variável de previsão em uma regressão, o modelo de médias-móveis usa erros de previsão do passado em um modelo de regressão-like.

$y_t = c + e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2} + \dots + \theta_q e_{t-q}.$

Esse modelo tem ordem $q$ e é denotado por $MA(q)$.

Para um modelo $MA(1)$, tem-se $\theta_1\in(-1,1)$

Para um modelo $MA(2)$, tem-se $\theta_1,\theta_2\in(-1,1)$ com $\theta_1 + \theta_2 > -1$ e $\theta_1 - \theta_2 < 1$

Relação entre modelos AR e MA

É possível escrever o modelo $AR(p)$ como um modelo $MA(\infty)$. Em particular para um modelo $AR(1)$, vem:

$y_t = \varphi_1 y_{t-1} + e_t$

$y_t = \varphi_1 ( \varphi_1 y_{t-2} + e_{t-1} ) + e_t$

$y_t = \varphi_1^2 y_{t-2} + \varphi_1 e_{t-1} + e_t$

$y_t = \varphi_1^3 y_{t-3} + \varphi_1^2 e_{t-2} + \varphi_1 e_{t-1} + e_t$

e assim por diante. Onde $\varphi_1\in(-1,1)$

Modelo ARIMA

ARIMA (Autoregressive Integrated Moving Average)

Combinação dos modelos:
- auto-regressivo (AR)
- médias-móveis (MA)
Método ARMA, acrescido da etapa de integração (I).
Caso a série temporal seja estacionária, o método se reduz ao ARMA.

Equação de Predição

\[y(t) = \sum_{i=1}^{p}\varphi_ix(t-i) + \sum_{j=1}^{q}\theta_j\epsilon(t-j)\]

Notação

$ARIMA(p,d,q)$, onde:

p = número de valores passados da variável (AR)
d = número de diferenças não-sazonais
q = número de rerros de predição passados (MA)

Método ARIMA em prospecções de médio prazo

Para 4 meses

mesesAPrever <- 4

Modelo:

arima_model_mensal <- auto.arima(tsMensalTrain)

Predição:

fcast_arima_4mth <- forecast(arima_model_mensal, h=mesesAPrever)
plot(fcast_arima_4mth, xlab="Anos", ylab="Fluxo Mensal", include=36)
lines(tsMensalTest4mth, col="red")

plot of chunk unnamed-chunk-41

Acurácia:

accuracy(fcast_arima_4mth, tsMensalTest4mth)

##                    ME  RMSE   MAE    MPE  MAPE   MASE      ACF1 Theil's U
## Training set    716.4 11391  8353  0.115 1.625 0.2089 -0.001948        NA
## Test set     -12192.3 18058 15355 -1.513 1.933 0.3840 -0.632429    0.2873

Para 6 meses

mesesAPrever <- 6

Predição:

fcast_arima_6mth <- forecast(arima_model_mensal, h=mesesAPrever)

plot(fcast_arima_6mth, xlab="Anos", ylab="Fluxo Mensal", include=36)
lines(tsMensalTest6mth, col="red")

Acurácia:

accuracy(fcast_arima_6mth, tsMensalTest6mth)

##                    ME  RMSE   MAE    MPE  MAPE   MASE      ACF1 Theil's U
## Training set    716.4 11391  8353  0.115 1.625 0.2089 -0.001948        NA
## Test set     -10206.2 15249 12315 -1.258 1.539 0.3080 -0.585487    0.2658

Método ARIMA em prospecções de longo prazo

Para 1 ano

mesesAPrever <- 12

Predição:

fcast_arima_1yr <- forecast(arima_model_mensal, h=mesesAPrever)

plot(fcast_arima_1yr, include=48)
lines(tsMensalTest1yr, col="red")

plot of chunk unnamed-chunk-49

Acurácia:

accuracy(fcast_arima_1yr, tsMensalTest1yr)

##                  ME  RMSE   MAE    MPE  MAPE   MASE      ACF1 Theil's U
## Training set  716.4 11391  8353 0.1150 1.625 0.2089 -0.001948        NA
## Test set     4503.2 18305 15764 0.4428 1.841 0.3943  0.403569    0.4769

Modelo TBATS

Técnica recente, publicada em 2011 (De Livera, Hyndman, and Snyder 2011).
A sigla TBATS é composta das iniciais das técnicas utilizadas: Trigonométrica, transformação Box-cox, correção de erro pela técnica ARMA, e componentes Trend e Sazonal.
Desenvolvido para séries temporais com características especiais de sazonalidade.

Justificativa

Por inspeção visual do fluxo diário, suspeitou-se que a série temporal apresenta sazonalidades semanal e anual.

Modelo TBATS, a partir da série multi-sazonal do conjunto de treinamento:

tbats_model <- tbats(tsDiarioTrain)
#tbats_model <- tbats(window(tsDiarioTrain, start=2005))

Decomposição do fluxo diário em termos trend+season pelo método TBATS

plot(tbats_model)

plot of chunk unnamed-chunk-52

O método TBATS é indicado para a previsão do fluxo diário, pois:

Sazonalidades múltiplas: semanal e anual
Sazonalidade semanal é de alta frequência

Método TBATS em prospecções de Curto Prazo

Para 30 dias

diasAPrever <- 30

Predição:

tbats_fc30Days <- forecast(tbats_model, h=diasAPrever)
plot(tbats_fc30Days, include = 120)
lines(tsDiarioTest30Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-54

Acurácia:

accuracy(tbats_fc30Days, tsDiarioTest30Days)

##                  ME  RMSE   MAE     MPE  MAPE   MASE       ACF1 Theil's U
## Training set -13.13 742.6 436.1 -0.1878 2.095 0.2186 -0.0007912        NA
## Test set     412.21 689.8 601.6  1.3682 2.122 0.3015  0.6355962    0.3331

Método TBATS em prospecções de Curto Prazo

Para 45 dias

diasAPrever <- 45

Predição:

tbats_fc45Days <- forecast(tbats_model, h=diasAPrever)
plot(tbats_fc45Days, include = 120)
lines(tsDiarioTest45Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-57

Acurácia:

accuracy(tbats_fc45Days, tsDiarioTest45Days)

##                  ME  RMSE   MAE     MPE  MAPE   MASE       ACF1 Theil's U
## Training set -13.13 742.6 436.1 -0.1878 2.095 0.2186 -0.0007912        NA
## Test set     572.57 869.5 732.4  1.8836 2.518 0.3670  0.6215786    0.3952

Exponential Smoothing

Médias móveis simples pondera igualmente as observações passadas em uma série temporal.
Exponential smoothing atribui pesos que decrescem exponencialmente ao longo do tempo, levando em consideração:
- Erro
- Tendência
- Sazonalidade
Estes parâmetros podem ser escolhidos como (e.g. tendência $T_h$, assumindo um termo de nível $l$ e termo de crescimento $b$, para $h$ períodos temporais e $0 < \phi < 1$):
- Nulo: $T_h = l$
- Aditivo: $T_h = l + bh$
- Aditivo amortecido: $T_h = l + \left( \phi + \phi^2 + \cdots + \phi^h \right)b$
- Multiplicativo: $T_h = l b^h$
- Multiplicativo amortecido: $T_h = l b^{\left( \phi + \phi^2 + \cdots + \phi^h \right)}$

Exponential Smoothing - 2

Ver (R. J. Hyndman and Athanasopoulos 2013). O autor possui um livro exclusivamente sobre o assunto em (R. Hyndman et al. 2008).

ETS

Para cada um dos 15 métodos, há dois modelos possíveis: com erro aditivo e com erro multiplicativo.

Há, portanto, 30 métodos distintos.

Exponential Smoothing no R

A função ets(ts) sem argumentos além da série temporal ts determina automaticamente o método apropriado.
O algoritmo avalia os 30 métodos e busca aquele que tiver o menor AIC.
- AIC - “Akaike information criterion”: uma estimativa da informação perdida quando um determinado modelo é usado para representar o verdadeiro processo que gera a série temporal.
Ao final, o objeto retornado por ets apresenta o modelo escolhido e o AIC resultante.
A função ets não trabalha com séries cuja sazonalidade é superior a 24 unidades temporais. Portanto, será utilizada somente pra prospecções realizadas com dados mensais (sazonalidade = 12).
O modelo ets parametrizado automáticamente deve dar resultados muito precisos para prospecções de poucos pontos (por exemplo, 4 pontos).

Método Exponential Smoothing em Médio Prazo

Para 6 meses

mesesAPrever <- 6

Modelo:

etsMensal <- ets(tsMensalTrain)

Prospecção:

fcastMensal6mth <- forecast(etsMensal, h=mesesAPrever)

Comparar predição com valores reais do conjunto de teste:

plot(fcastMensal6mth, xlab="Anos", ylab="Fluxo Mensal", include = 36);
lines(tsMensalTest6mth, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-62

Acurácia:

accuracy(fcastMensal6mth, tsMensalTest6mth)

##                  ME  RMSE  MAE      MPE  MAPE   MASE     ACF1 Theil's U
## Training set  856.7  9849 7308  0.13394 1.442 0.1828 -0.02704        NA
## Test set     -755.7 10937 8653 -0.06425 1.104 0.2164 -0.63165    0.2848

Prospecção a Longo Prazo

Transformação Box-Cox

Há um fator multiplicativo no fator sazonal. Uma transformação Box-Cox demonstrou ser apropriada.

$W_t = \begin{cases} \log\left(y_t\right) & \lambda = 0\\ \frac{(y_t^\lambda - 1)}{\lambda} & \lambda \neq 0\end{cases}$

Série original

plot(tsMensalTrain, xlab="Anos", ylab="Fluxo Mensal")

plot of chunk unnamed-chunk-64

Série transformada

Calcula-se o fator $\lambda$ para a transformação:

lam <- BoxCox.lambda(tsMensalTrain)

Série Resultante:

tsMensalBoxCox <- BoxCox(x = tsMensalTrain, lam)

plot(tsMensalBoxCox, col="red")

plot of chunk unnamed-chunk-67

Prospecção com decomposição sazonal de Loess - STLF

Aplica-se a decomposição STL
Remove-se a sazonalidade da série temporal
Um modelo é treinado na série resultante
Faz-se a prospecção
Adiciona-se novamente o último período da sazonalidade estimada nos resultados

mesesAPrever <- 24

Dados com sazonalidade removida:

tsMensal.stl <- stl(tsMensalTrain[,1], s.window=12)
# Seasonally adjusted data constructed by removing the seasonal component.
plot(seasadj(tsMensal.stl))

plot of chunk unnamed-chunk-69

stlf_model <- stlf(tsMensalTrain[,1], lambda = lam)
stlf_fcast <- forecast(stlf_model, h=mesesAPrever)
plot(stlf_fcast, xlab="Anos", ylab="Fluxo Mensal")
lines(tsMensalTest2yr, col="red")

plot of chunk unnamed-chunk-70

Acurácia:

accuracy(stlf_fcast, tsMensalTest1yr)

##                   ME  RMSE   MAE      MPE  MAPE   MASE     ACF1 Theil's U
## Training set  -25.42  8186  6175 -0.01584 1.238 0.1544 -0.02348        NA
## Test set     8132.24 17138 13696  0.92760 1.619 0.3425 -0.11417    0.4831

Referências

De Livera, Alysha M, Rob J Hyndman, and Ralph D Snyder. 2011. “Forecasting Time Series with Complex Seasonal Patterns Using Exponential Smoothing.” Journal of the American Statistical Association 106 (496): 1513–1527.

Hyndman, Rob J. 2014. “Forecasting Time Series Using R.” http://robjhyndman.com/talks/MelbourneRUG.pdf.

Hyndman, Rob J, and George Athanasopoulos. 2013. Forecasting: principles and Practice. OTexts.

Hyndman, Rob, Anne B. Koehler, J. Keith Ord, and Ralph D. Snyder. 2008. Forecasting with Exponential Smoothing: The State Space Approach (Springer Series in Statistics). 2008th ed. Springer. http://amazon.com/o/ASIN/3540719164/.

Leek, Jeffrey. 2014. “Practical Machine Learning - Lecture 27 - Forecasting.” Johns Hopkins Bloomberg School of Public Health. https://d396qusza40orc.cloudfront.net/predmachlearn/027forecasting.pdf.

Prospecção Séries Temporais

Importação da base de dados

Séries Temporais

Série Mensal:

Série Diária:

Inspeção dos Dados

Série para um dia específico da semana

Série para um dia e mês específico da semana

Inspeção dos Dados - 2

Fluxo Mensal

Fluxo Diário

Observações

Inspeção dos Dados - 3

Fluxo Diário para um intervalo arbitrário de 8 semanas

Inspeção dos Dados - 4

Auto-correlação

Parcial:

Objetivos

Curto Prazo

Médio Prazo

Longo Prazo

Separação dos dados

Curto Prazo (tsDiario)

Conjunto de Treinamento

Conjuntos de Teste

30 dias:

45 dias:

Médio/Longo Prazo (tsMensal)

Conjunto de Treinamento

Conjuntos de Teste

4 meses:

6 meses:

1 ano:

2 anos:

Bibliotecas Utilizada

Prospecções Estatísticas

Redes Neurais:

Estratégia de Avaliação dos Métodos

Acurácia

MAPE:

Benchmark Métodos Estatísticos

Curto prazo

Médio prazo (4 meses)

Longo prazo

Modelo auto-regressivo (AR)

Modelos médias-móveis (MA)

Relação entre modelos AR e MA

Modelo ARIMA

Equação de Predição

Notação

Método ARIMA em prospecções de médio prazo

Para 4 meses

Para 6 meses

Método ARIMA em prospecções de longo prazo

Para 1 ano

Modelo TBATS

Justificativa

Modelo TBATS, a partir da série multi-sazonal do conjunto de treinamento:

Decomposição do fluxo diário em termos trend+season pelo método TBATS

Método TBATS em prospecções de Curto Prazo

Para 30 dias

Método TBATS em prospecções de Curto Prazo

Para 45 dias

Exponential Smoothing

Exponential Smoothing - 2

Há, portanto, 30 métodos distintos.

Exponential Smoothing no R

Método Exponential Smoothing em Médio Prazo

Para 6 meses

Prospecção a Longo Prazo

Transformação Box-Cox

Série original

Série transformada

Prospecção com decomposição sazonal de Loess - STLF

Referências

Curto Prazo (`tsDiario`)

Médio/Longo Prazo (`tsMensal`)